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ABSTRACT 

In this paper we present a critical discussion of the algorithms commonly used in 
N-body simulations of Galaxy Formation to deal with the energy equation governing 
heating and cooling, to model star formation and the star formation rate, and to 
account for the energy feed-back from stars. 

First, we propose our technique for solving the energy equation in presence of 
heating and cooling, which includes some difference with respect to the standard semi- 
implicit techniques. 

Second, we examine the current criteria for the onset of the star formation activity. 
We suggest a new approach, in which star formation is let depend on the total mass 
density - baryonic (gas and stars) and dark matter - of the system and on the metal- 
dependent cooling efficiency 

Third, we check and discuss the separated effects of energy (and mass) feed-back 
from several sources - namely supernovae, stellar winds from massive stars, and UV 
flux from the same objects. 

All the simulations are performed in the framework of the formation and evolution of 
a disk galaxy. 

We show that the inclusion of these physical phenomena has a significant impact on 
the evolution of the galaxy models. 

Key words: Numerical methods: SPH; Galaxy: formation, evolution; Stars: forma- 
tion, feed-back 



1 INTRODUCTION 

Numerical hydrodynamics and N-body simulations are 
nowadays the fundamental tool to investigate galaxy for- 
mation and evolution. The current status of numerical and 
semi-analytical modeling of various galaxy properties has 
been recently summarized by Frenk et al. (1997). 

One of the poorly understood processes is the forma- 
tion of stars. Since the pioneering study by Katz (1992), 
the onset of star formation (SF) is usually empirically pa- 
rameterized (see also Gerritsen 1997). The key idea is that 
an element of fluid must satisfy some conditions in order to 
be eligible to SF and turn part of its gas content into stars 
according to a suitable star formation rate (SFR), which is 
customarily a reminiscence of the Schmidt (1959) law. The 
reader is referred to the exhaustive discussion of this topic 
by Mihos & Hernquist (1994) for more details. 

The basic criteria to select the fluid elements prone to 
SF are (i) the gas particle must be in a convergent flow, and 
(ii) the gas particle must be Jeans unstable. The conditions 



are met if the velocity divergence is negative and if the sound 
crossing time-scale is shorter than the dynamical time-scale. 

Schmidt-like laws imply that SFR oc p n , where the 
power n ranges between 1 and 2. Conservation arguments 
suggest that the most probable value of n is § (see Katz 
1992). The Schmidt law written for the volume density is 



dp-,, 
~~dt 



dpg 
dt 



c*£g 

to 



(1) 



where is the so-called dimensionless efficiency of star for- 
mation, and t g is the characteristic time for the gas to flow. 
This is chosen to be the maximum between the cooling time 
and the free-fall time tff = (4nGp g )~z . 

Over the past few years, starting from this general back- 
ground a number of refinements have been brought forth. 
These mainly concern the criteria to decide whether a gas- 
particle is suitable to form stars, and the way stars emerge 
from gas, i.e. how gas clouds undergo fragmentation. 

Navarro & White (1993) and Katz, Weinberg & Hern- 
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quist (1996) introduce an additional condition for the onset 
of SF, i.e. the so-called over-density criterion, which secures 
that a collapsing cloud remains cool. The threshold density 
implies that the cooling time-scale is much shorter than the 
dynamical time-scale, and this always occurs for tempera- 
tures greater than the cut-off temperature of the cooling 
curve (T = W 4 K for the typical mass resolution of this type 
of problem). In addition to this, Katz, Weinberg & Hernquist 
(1996) impose that the gas-particles can form stars only if 
a minimum physical density corresponding to 0.1 hydrogen 
atoms per cm 3 is met. 

In Katz et al. (1996) and Navarro & White (1993), the 
fragmentation of a gas cloud is usually modeled in the follow- 
ing way. When in a gas-particle SF switches on, the gaseous 
mass decreases until it reaches some minimum value, say 
5% of the original value. At this stage, the gas-particle is 
turned into a collision-less star-particle, while the remain- 
ing gas is distributed among the surrounding gas-particles. 
In the meantime the gas-particle is considered as a dual en- 
tity, whose mass is partly collision-less and partly collisional. 
This approach is a good compromise in terms of computa- 
tional effort because it is not necessary to add a collision-less 
particle every time a new star is formed (Katz, Weinberg & 
Hernquist 1996). Then the choice of the star forming regions 
is made by means of suitable probability arguments. Every 
time a gas-particle forms stars, its mass is reduced by 1/3 
(free parameter) . A different value is adopted by Navarro & 
White (1993), who suppose that at each star forming event 
the gas mass is halved, and that each gas-particle can split 
into four star-particles at most (this limit is set by compu- 
tational limitations). 

A somewhat different scheme is followed by Steinmetz 
(1995) who calculates the mass of the newly born star- 
particle using equation (1) in the form: 

m* = m s (l - exp(-^^)), (2) 
iff 

where At is the particle time step. The mass of the gas- 
particle is accordingly reduced. Although this latter ap- 
proach appears to be more physically sounded, compu- 
tational problems can arise when many star-particles are 
expected to form, so that the number of star formation 
episodes must be artificially limited. 

Once stars are present, they are expected to return to 
the interstellar medium (ISM) part of their mass (in form of 
chemically processed gas) and energy via Supernovae (SNae) 
explosions, stellar winds, and UV flux. These latter contri- 
butions are significant only in the case of massive stars. All 
this is known as the stellar energy Feed-Back (FB). 

The amounts of mass, metals and energy released by 
each star depend on a star's mass, lifetime, and initial com- 
position in a way that can be easily taken into account. 

The key problem here is to know how the energy re- 
leased by stars is given to the ISM because the limited reso- 
lution of N-body simulations does not allow one to describe 
the ISM as a multi-phase medium. Basically two different 
schemes exist: 

(i) all the energy (from SNae, stellar winds, and UV 
flux) is given to the thermal budget of the fluid element 
(Katz 1992, Steinmetz & Miiller 1994); 

(ii) only a fraction of this energy is let kinetically affect 
the surrounding fluid elements. This fraction (/„ or ekin) is a 



free parameter (Navarro & White 1993, Mihos & Hernquist 
1994) chosen in such a way that some macroscopic properties 
of the galaxy, like for instance the metallicity distribution 
(Groom 1997) or the shape of the disk (Mihos & Hernquist 
1994), are reproduced. 

The deposit of all the energy in the thermal budget 
of the fluid (first alternative) has little effect on a galaxy' 
structure and evolution. In fact, since the thermal energy is 
given to a medium of high density and short cooling time- 
scale in turn, it is almost immediately radiated away. The 
difficulty cannot be cured by supposing that the thermal 
energy is injected and shared over an e-folding time scale 
(Summers 1993). 

In all the above schemes the fraction of energy from 
stellar FB which is not transformed into kinetic energy is 
added to the energy budget of the fluid element, and the 
fluid element is let cool. The integration of the energy equa- 
tion with respect of time is generally made according to 
a semi-implicit scheme (see the discussion in Hernquist & 
Katz 1989). 

In this paper we present in detail our prescriptions for 
computing energy balance, SF and FB. The organization of 
the paper is as follows. Section 2 describes the numerical 
code and the method used to integrate the energy equation, 
whereas Section 3 discusses in details the adopted initial con- 
ditions. In section 4 we test our energy integration scheme 
following the formation of a galactic disk. Section 5 presents 
tests a discussion on various SF criteria , our new approach 
to SF and the corresponding models. Section 6 examines 
the effects of different sources of FB, while, finally, some 
concluding remarks are presented in Section 7. 

2 THE BASIC NUMERICAL TOOL 
2.1 The Tree-SPH code 

All the simulations presented here have been performed us- 
ing the Tree-SPH code developed by Carraro (1996) and Lia 
(1996), and described by Carraro et al. (1998a), to whom the 
reader is referred for all details. 

The code, which can follow the evolution of a mix of 
Dark Matter (DM) and Baryons (gas and stars), has been 
carefully checked against several classical tests with satis- 
factory results, as reported in Carraro et al. (1998a). 

In this code, the properties of the gas component are 
followed by means of the Smoothed Particle Hydrodynam- 
ics (SPH) technique (Lucy 1977, Gingold & Monaghan 1997, 
Benz 1990), whereas the gravitational forces are computed 
by means of the hierarchical tree algorithm of Barnes & Hut 
(1986) using a typical tolerance parameter 8 = 0.8 and ex- 
panding tree nodes to quadrupole order. We adopt a Plum- 
mer softening parameter. 

In SPH each particle represents a fluid element whose 
position, velocity, energy, density etc. are followed in time 
and space. The properties of the fluid are locally estimated 
by an interpolation which involves the smoothing length hi. 
In our code each particle possesses its own time and space 
variable smoothing length hi , and evolves with its own time- 
step. This renders the code highly adaptive and flexible, and 
suited to speed-up the numerical calculations. 

Radiative cooling is described by means of numerical 
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tabulations as a function of temperature and metallicity 
taken from Sutherland & Dopita (1994). This allows us to 
account for the effects of variations in the metallicity among 
the fluid elements and for each of these as a function of time 
and position. 

The chemical enrichment of the gas-particles caused by 
SF and stellar ejecta (Portinari et al 1998) is described by 
means of the closed-box model applied to each gas-particle 
(cf. Carraro et al. 1998a,b for more details). 

Star formation and Feed-back are discussed below. 

Finally, all the calculations presented here have been 
carried out on a DIGITAL ALPHA-2000 (330 Mhz of Clock, 
512 Mbyte of RAM) workstation hosted by the Padua Ob- 
servatory & Astronomy Department. 



2.2 On the integration of the energy equation 

The usual form of the energy equation in SPH formalism is 

\ {ViWimM)* VtWinj,^)) + r ~ Ao , (3) 

2 v ' p 

(Benz 1990; Hernquist & Katz 1989). The first term rep- 
resents the heating or cooling rate of mechanical nature, 
whereas the second term F is the total heating rate from all 
sources apart from the mechanical ones, and the third term 
Ac I p is the total cooling rate by many physical agents (see 
Carraro et al. 1998a for details). 

In absence of explicit sources or sinks of energy the 
energy equation is adequately integrated using an explicit 
scheme and the Courant condition for time-stepping (Hern- 
quist & Katz 1989). 

The situation is much more complicated when consider- 
ing cooling. In fact, in real situations the cooling time-scale 
becomes much shorter than any other relevant time-scale 
(Katz & Gunn 1991), and the time-step becomes consid- 
erably shorter than the Courant time-step, even using the 
fastest computers at disposal. This fact makes it impossible 
to integrate the complete system of equations (cf. Carraro 
et al. 1998a) adopting as time-step the cooling time-scale. 

To cope with this difficulty, Katz & Gunn (1991) damp 
the cooling rate to avoid too short time-steps allowing gas 
particles to loose only half of their thermal energy per time- 
step. 

Hernquist & Katz (1989) and Dave et al. 1997 solve 
semi-implicitly eq. (||) using the trapezoidal rule, 

= Ui + -{dti x e; + dt^ xe^), (4) 

where = dui/dt. The leap-frog scheme is used to update 
thermal energy, and the energy equation, which is nonlin- 
ear for Uj, is solved iteratively both at the predictor and 
at the corrector phase. The technique adopted is a hybrid 
scheme which is a combination of the bi-section and Newton- 
Raphson methods (Press et al 1989). The only assumption 
is that at the predictor stage, when the predicted u" +1 is 
searched for, the terms u" +1 are equal to it™ -1 . 



Our scheme to update energy is conceptually the same, 
but differs in the predictor stage and in the iteration scheme 
adopted to solve equation ^[ 

In brief, at the first time-step the quantity u n ~ 1 ' 2 is 
calculated and for all subsequent time steps, the leap-frog 
technique, as in Steinmetz & Miiller (1993), is used: 

(i) We start with u™ at t n ; 

(ii) compute u n as 

~n n — 1/2 . 1 , n , n 

u = u ' + — t x e . 

This predicted energy, together with the predicted velocity 
is used to evaluated the viscous and adiabatic contribution 
to e™ . In other words the predictor phase is calculated 
explicitly because all the necessary quantities are available 
from the previous time step t n . 

(iii) finally, derive u n+1 solving the equation ^ itera- 
tively (corrector phase) for both the predicted and old adi- 
abatic and viscous terms; 

In the corrector stage the integration of the equation ^ is 
performed using the Brent method (Press et al 1989) instead 
of the Newton-Raphson, the accuracy being fixed to a part 
in 10 . The Brent method has been adopted because it is 
better suited as root-finder for functions in tabular form 
(Press et al. 1989). 



3 INITIAL CONDITIONS 

In the following we are going to check our treatment of the 
energy equation and our implementation of SF and FB sim- 
ulating the formation of a disk galaxy. Since our code is not 
a cosmological code we shall consider an ad hoc initial con- 
figuration for a protogalaxy. This approach is justified by 
the fact they we are going to focus on processes occurring 
at scales much lower than the cosmological ones. 



3.1 Theoretical overview 

According to modern gravitational instability theory (White 
& Rees 1978), protogalaxies comprise a mixture of 
dissipation-less dark matter (DM), whose nature is not clear 
insofar, and dissipational baryonic material, roughly in the 
mass ratio 1:10. After a violent relaxation process, which fol- 
lows the separation from the Hubble expansion, a DM halo 
becomes isothermal, and acquires baryonic material which 
heats up at the halo virial temperature. Gas then cools and 
collapses, and through fragmentation turns into the stars we 
see today in disks and spheroidal systems. 

In this scenario disks form as a consequence of the an- 
gular momentum that DM halos acquire due to the tidal 
torque felt by surrounding halos, while spheroidal systems 
are thought to be produced by merging of disks. 

Instead of selecting halos from cosmological N-body 
simulations (see for instance Weil et al 1998) , in our code we 
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Table 1. Initial conditions of the simulated Dark Matter halo: 
N=10,000: a:b:c = 3.0: 2.25: 1.5. 
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Figure 1. The initial specific angular momentum j as a function 
of the radial distance r from the center of a 1/r dark matter halo. 



set up a protogalaxy as an isolated rotating DM halo with 
baryonic material inside, proceeding as follows. 



3.2 Initial configuration 

We consider a triaxial DM halo whose density profile is 

p(r) oc -. (5) 
r 

Although rather arbitrary, this choice seems to be quite 
reasonable. Indeed DM halos emerging from cosmological 
N-body simulations are not King or isothermal spheres, 
but show, independently from cosmological models, initial 
fluctuations spectra and total mass, an universal profile 
(Navarro et al 1996; Huss et al 1998). This profile is not 
a power law, but has a slope a = dlnp/dlnr with a = — 1 
close to the halo center, and a — —3 at larger radii. Thus in 
the inner part the adopted profile matches the universal one. 
Moreover this profile describes a situation which is reminis- 
cent of a collapse within an expanding universe, being the 
local free fall time a function of the radius (see the discus- 
sion in Curir et al 1993; Aguilar & Merritt 1990 ). Moreover 
triaxiality is quite natural for galactic DM halos (Becquaert 
& Combes 1997; Oiling & Merrifield 1998). 

DM particles are distributed by means of an acceptance- 
rejection criterium, following step by step the procedure de- 
veloped by Curir et al (1993). Specifically, we consider a 
triaxial ellipsoid whose radial scale r is given by 



A 


P 


Si 


<5 2 


0.09 


0.20 


0.05 


0.025 



where a > b > c are the axial ratios (see Table 1). 
The velocity field has been chosen to produce an angular 
momentum which depends linearly on the distance, j(r) oc r 
(Barnes & Efstathiou 1987). It has been built up sampling 
the moduli V of the particles velocity by using a Maxwellian 
distribution and imposing that the virial ratio j3 — 2xT/\W\ 
is equal to some fixed value ranging between 0.05 and 0.20. 
Here T is the kinetic energy, while W is the potential energy. 
The final cartesian components of the particle velocity can 
be derived from the following formulas 



V x = -V x sin(9 + a) 
V v = Vx cos(6 + a) 



(7) 
(8) 
(9) 



(ac + 



V 



+ 



(b/a)> (c/a 



(6) 



where 9 is the angle between the position vector of a parti- 
cle and the x coordinate axis, a is an angle varying between 
— Sin and +<5i7r, whereas £ is a random parameter ranging 
from —82V and +82V . 81 and 82 are chosen so that the ki- 
netic energy T changes by less than 1% of the value giving 
the initial virial ratio f3. 

Table 1 summarizes the adopted initial values for the set of 
parameters introduced above. 

Due to the assigned velocity field, the halo acquires an 
amount of angular momentum, which is conventionally de- 
scribed by means of the dimensionless spin parameter A: 

x _ J\E\ 1/2 
CM 5 /? 

Here G is the gravitational constant, J the system angular 
momentum, M the total mass and E the total system energy. 
In our case the A parameter has chosen to be 0.09, signif- 
icantly greater than the mean values of cosmological halos 
(Steinmetz & Bartelmann 1995). This choice is motivated 
by the comparison we are going to make with similar initial 
conditions (Navarro & White 1993; Thacker et al 1998; Rai- 
teri et al 1996). 

The softening parameter e is computed as follows. Af- 
ter plotting the inter-particles separation as a function of 
the distance to the model center, we compute e as the mean 
inter-particles separation at the center of the sphere, tak- 
ing care to have at least one hundred particles inside the 
softening radius (Romeo A. C, Pearce F. R., private com- 
munications). We consider a Plummer softening parameter, 
equal for both DM and gas particles (see below), and keep 
it constant along the simulation. For this particular choice 
of the initial configuration the softening parameter e turns 
out to be 3.6 kpc. 

Total energy and angular momentum are conserved 
within 1% and 0.1% level, respectively. 

We consider a dark matter halo with mass 1O 12 M and 
radius (a) 120 kpc in order to simulate a galaxy with a size 
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Figure 2. The formation process of a galactic disk in the x — z plane. In the bottom— right corner of any snapshot time is reported in 
billion years. 



similar to the Milky Way. Axial ratios are a = 3.00, b = 2.25, 
c= 1.50. 



4 THE FORMATION OF A GASEOUS DISK 

The formation of a galactic disk is simulated distributing 
gas inside the halo described above and switching cooling 
on. To mimic the in-fall of gas inside the potential well of 
the halo, we distribute gas particles (10,000 in number) on 
the top of DM particles. The baryonic fraction adopted is 
fb — 0.1 , and gas particles are Plummer-softened in the 
same way as the DM particles. 

Under cooling and the velocity field of the halo, the gas 
is expected to settle down in a rotating thin structure. For 



this purpose we need to specify for the baryonic component 
a temperature and a metal content. We assume that the gas 
has a temperature of « 1.0 x 10 4 °K (Navarro & White 
1993), and an almost primordial metal content amounting 
to Z ~ 10~ 4 , which translates into [Fe/H] ~ — 3 (Bertelli 
et al 1994). 

This simulation is meant to verify that our energy inte- 
gration scheme works properly when cooling is switched on. 
For this reason SF, FB and chemical enrichment are turned 
off. 

The formation of the disk is shown in Fig. 3. Left panels 
follows from the top to the bottom the evolution of the gas 
component in the x — z plane. Time is in Gyr. 
Starting from an ellipsoidal initial configuration (t = 0.0), 
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Figure 3. The Star Formation history ( in solar masses per year) for four models in which star formation criteria are defined as in 
Section. 5. See text for details 



baryons due to cooling and angular momentum settle down 
in a thin rotating structure (t = 1.75), eventually developing 
a bulge like structure in the model center (t = 2.5). At this 
time the gaseous disk is about 5 kpc high and 40 kpc wide. 
Comparing our results with similar simulation by Raiteri et 
al (1996) we see that in our case the disk forms more slowly. 
This is due to the lower efficiency of our cooling with respect 
to the more widely adopted Katz & Gunn (1991) cooling. 



5 STAR FORMATION 

As already recalled, SF in N-body simulations is let occur 
if suitable criteria meant to simulate the real behaviour of 
the ISM are satisfied. 

Current prescriptions for SF stand on three time-scales, 



i.e. the crossing, cooling, and free-fall time-scales, which ul- 
timately depend (although in a different fashion) on the den- 
sity of the fluid. According to the most popular prescription 
(Katz 1992; Navarro & White 1993), SF is let occur if 

* V ■ Vi < 

* tsound ^* "ff 

* tcooling "ff 

Aim of this section is to closely scrutiny the effects of three 
conditions on the model results. We start with some general 
considerations. 

Most likely, the divergence criterion is not strictly re- 
quired. Firstly the typical mass resolution (10 6 - 1O 7 M ) 
of this kind of N-body simulation is clearly not sufficient to 
follow the kinematical behaviour of the gas at the molecular 
clouds scale. Then in principle a fluid element can form stars 
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Table 2. Comparison of SF criteria. 



Case 


SF peak 


Peak time 


Stars formed 




M Q /yr 


Gyr 


Number 


A 


85 


2.25 


3381 


B 


73 


2.10 


3771 


C 


71 


2.03 


6580 


D 


77 


1.81 


8543 



without being in a convergent flow (let us think about shock 
induced or stochastic SF ). For instance Mihos & Hernquist 
(1994) adopt in the special case of an already formed spi- 
ral galaxy a different dynamical criterium, based upon the 
Toomre instability criterium. 

The free-fall time scale is usually computed as 

tff = (4nGp)- 1/2 (10) 

where p is the gas density of the fluid element. We have 
replaced this condition with the more general one 

t ff = (4irG Ptot )- 1/2 (11) 

where p to t = pgas + Pdm- In principle, one should consider 
also p s tar- In fact the dynamical behaviour of the fluid el- 
ement under consideration depends on all the mass inside 
the volume with radius 2 x h, where h is the smoothing 
length. Therefore at the typical resolution of these simula- 
tions (about 2.5 kpc) the presence of DM accelerates the 
collapse of a fluid element by lowering the free-fall time- 
scale. DM density is evaluated in the same manner as gas 
density, defining a smoothing length for any DM particle, 
although this scheme results computationally expensive. 
The importance of DM on SF has been recently studied from 
an analytical point of view by Caimmi & Secco (1997). 

Finally, we remind the reader that the condition on the 
cooling time-scale is customarily replaced by a condition on 
the density, i.e. a constant density threshold (Navarro & 
White 1993). This stems from suitable arguments about the 
cooling time-scale based on the notion that the gas metallic- 
ity has no effect. However, as the cooling rate does increase 
with the metal content of the ISM with consequent increase 
of the threshold density, the more general condition on the 
cooling time-scale ought to be preferred. 



5.1 Testing the SF recipes 

To this aim we perform the same simulation as in Section 4, 
say the formation of a disk galaxy, but letting stars form 
according to different set of criteria. We follow the model 
till the end to the major star formation episode. 
The simulations have been made using 10,000 DM and 
10,000 baryonic particles, neglecting in this particular set 
of models any source of FB. The following four cases are 
examined. 
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Figure 4. The heating rates from the three sources of FB as 
indicated. See the text for more details. 



In our models, fragmentation of gas-particles undergoing SF 
is let occur and new star-particles are created when the gas 
content of the fluid element has fallen below 20% of the 
initial value, being the remaininig gas spread out over the 
surrounding particles. 

This limits the number of star particles produced. More- 
over the gas particle is allowed to experiment up to ten SF 
episodes. Afterwards it cannot make stars anymore. 
In all simulations the efficiency parameter c* is fixed to 
0.1. Finally, when more than two star-particles form in- 
side a sphere with radius equal to the softening parameter 
e = 3.6 kpc, they are merged together to form a single ob- 
ject, in order to keep the total number of particles small. 

The results for the SFR as a function of time are shown 
in Fig. 3, whereas a few relevant quantities of the models 
are given in Table 2. 

Comparing either case A to case B or case C to case D, we 
notice that the criterion on the velocity divergence has in 
practice no effect, as argued above. 

In contrast, there is a significant difference passing from the 
cases A and B based on the over-density condition to the 
cases C and D based on the cooling time-scale. In fact the 
models of the first group start sensibly later to form stars 
and after the burst, which occurs somewhat later, SF main- 
tains higher that in the second group models, and show some 
secondary peaks. Dropping the divergence criterium signifi- 
cantly increases the number of spawned stars (model C and 
D). 

In these models after the peak SF proceeds with a rate 
around 4 — 6 MQ/yr, quite usual for disk galaxies. Putting 
together all these differences and the discussion above we 
are led to prefer the combination of criteria of the model 
D, which we are going to use in all the other simulations 
presented in this paper. 
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All simulations exhibit the SF peak at about 2 Gyr, much 
later than, for instance, the Raiteri et al (1996) simulation. 
This is clearly due to the different cooling functions adopted. 
Katz & Gunn (1991) analytical formulas provide a much 
more efficient cooling (see below and Carraro et al 1998a). 



6 STELLAR FEED-BACK 

In this section, we test separately the effect of FB of differ- 
ent nature making use of the case D for the recipe of star 
formation. Furthermore, we include the effect of chemical 
enrichment. 

The energy FB originates from SNae explosions of both 
type la and II (Greggio & Renzini 1983), the stellar winds 
from massive stars (Chiosi & Maeder 1986) and the UV flux 
from the same objects (Chiosi et al. 1998). 

In order to evaluate the amount of energy injected 
into ISM by the above sources we suppose that each newly 
formed star-particle, which for the mass resolution of our 
simulations has a total mass of the order of 10 7 M©, is actu- 
ally made of a larger number of smaller subunits (the real 
stars) lumped together and distributed in mass according to 
a given initial mass function (IMF). For the purposes of the 
present study we adopt the IMF by Miller & Scalo (1979) 
over the mass interval from 0.1 to 120 M©. 

Chemical evolution is followed by means of the closed- 
box approximation (Tinsley 1980; Portinari et al 1998) and 
sharing of the metals among the gas-particles is described 
by means of a diffusive scheme (Groom 1997; Carraro et al. 
1998a). 

Details on implementation of FB in our Tree-SPH code 
has been already reported in Carraro et al. (1998a), to whom 
the reader should refer. With respect to the previous study, 
only two major changes have been made. First, of the UV 
flux emitted by massive stars only a small fraction (0.01) 
is actually used. All the remaining UV flux is supposed to 
be re-processed by dust into the far-infrared and lost by 
the galaxy (Silva et al 1998). Second, the energy released to 
the ISM by a SN event amounts only to 10 50 ergs, roughly 10 
times less than in older evaluations (see Thornton et al. 1998 
for an exhaustive discussion of this topic) . The heating rates 
from the three sources of energy are displayed in Fig. 4 for 
the sake of illustration. A notable feature to remark, is that 
the energy injection by stellar winds may parallel that from 
SNae explosions, and the contribution from the UV flux plays 
a significant role although a strong dumping factor has been 
applied. All the FB mechanisms have almost the same trend, 
since they derive basically from the same stellar sources. 

Finally, we stress that all the energy from FB is given to 
the thermal budget of the particles, because the resolution 
we are working with makes it impossible to describe this 
process in more detail. In fact, the space resolution of these 
simulations is about 2.5 kpc, much larger than the typical 
distance over which the effects of SNae explosions and stel- 
lar winds are visible (Mckee & Ostriker 1977). See also the 
arguments given by Carraro et al. (1998a). Possible interac- 
tions of kinetic nature among gas particles occur only via 
the pressure gradients. 



6.1 Testing the sources of Feed-Back 

In order to understand the role played by each source of FB 
we perform the following experiments using the same model 
as in the previous section. Three cases are considered: 

Case Dl: FB from SNae; 

Case D2: FB from SNae and stellar winds; 

Case D3: FB from SNae, stellar winds, and UV flux. 

First of all in Fig. 5, for the case Dl, we show the SF history, 
and the rate of SNae of both types, la and II. Due to the dif- 
ferent stellar progenitors, SNae la and II rates show different 
trends. SNae II explode immediately, and their rate strictly 
follow the SF history. On the contrary, SNas la start to ap- 
pear later, being the progenitor less massive. The relative 
abundance of SNae events clearly depends on the adopted 
IMF. 

The results of these simulations are shown in the vari- 
ous panels of Fig. 6. A few relevant quantities of the models 
are given in Table 3. Namely we consider the mean gas tem- 
perature, the SF peak and the time at which it occurs, the 
final maximum metal abundance, and the radial and vertical 
dimensions of the final stellar disk. 

Looking at these results the following conclusions can be 
drawn. The models Dl and D2 do not differ significantly, 
simply turning on FB from stellar winds the final models 
becomes a bit hotter, less stars are formed and the final 
maximum metal abundance achieved is slightly lower. Since 
almost the same amount of stars are generated, the final 
stellar disk morphology is basically identical. 
When the combined effect of all the FB sources are consid- 
ered (model D3) the situation changes more significantly 
(see also Fig. 4). The model becomes hotter more rapidly 
when SF starts, SF peak occurs before, and less stars arc 
formed. As a consequence the final mean temperature is 
somewhat lower and the disk morphology is not well resolved 
as in the previous models. The final maximum metallicity 
turns out somewhat greater due to the effect of the diffu- 
sive mixing. Metals are spread and smooth around the gas 
neighbors. Since SF in this model stops before, no metal 
poor stars are produced, and the mean metallicity does in- 
crease. 

In all the models a number of gas particles do not ex- 
perience SF and keep the primordial metal content. 

7 SUMMARY AND CONCLUSIONS 

In the framework of the formation of a spiral-like galaxy, we 
have examined in detail the effect of different prescriptions 
for Star Formation and Feed-Back on numerical simulations 
of galaxy formation and evolution. 

First of all we described in details how we integrate en- 
ergy equation when cooling is switched on. We adopted cool- 
ing functions which depend on metal abundance, and use 
the Brent method as root finder instead of the more widely 
used Newton-Raphson scheme. Brent scheme performs bet- 
ter when using functions in tabular form. 

As far as the physical conditions at which Star Forma- 
tion is likely to start, we argue that the criterion on the 
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Figure 5. SF history and SNae rates for the model D3. SNas II closely follow SF history, whereas SNae la exhibit a peak somewhat later, 
as expected from the different progenitors lifetime. 



Table 3. Results for models with different Feed-Back. 



Case 


Tave 


SF peak 


Peak time 


[Fe/ H]f in 


rn 


I'D 




10 6 °K 


M e /yr 


Gyr 


dex 


kpc 


kpc 


Dl 


1.15 


39 


1.72 


0.446 


30 


5 


D2 


1.46 


38 


1.72 


0.406 


30 


5 


D3 


1.33 


32 


1.67 


0.513 


20 


5 



velocity divergence has no sizable effect on the model re- 
sults. 

Moreover the over-density and cooling time criteria are 
equivalent only if the effect of metallicity on the cooling 
time-scale is neglected. However, as in real galaxies the 
metallicity varies as a function of time and space, the crite- 



rion based on the cooling time-scale ought to be preferred. 
This finding is a step further toward better understanding 
under which conditions Star Formation takes place in real 
galaxies. 

In addition we introduce in the gas free-fall time computa- 
tion the contribution of DM and eventually stars which lie 
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Figure 6. Comparison of the galactic models evolution by using different sources of Feed-Back. The top row shows the evolution of a 
model in which only SNa3 FB is considered, middle row a model in which also stellar winds are switched on, and bottom row the case of 
a model in which all the FB sources are allowed to act. Right panels show the chemical enrichment (solid line the maximum metallicity, 
dotted line the minimum), middle panels the SF history, and left panels the final morphology of the stellar disk. 



within the neighbors sphere. This has the effect to decrease 
the free-fall time, accelerating the cloud collapse. 

Finally we have quantitatively shown the separated and 
cumulated effects of different sources of Feed-Back (SNag^ 
stellar wind and UV flux from massive stars) on the global 
properties (SFR, metallicity and morphology) of galaxy 
models. 
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